/*
 * EquilibriumSolver.h
 *
 *  Created on: 26 Jul 2011
 *      Author: Allan
 */

#ifndef EQUILIBRIUMSOLVER_H_
#define EQUILIBRIUMSOLVER_H_

// C++ includes
#include <vector>
#include <functional>
using namespace std;

// Eigen includes
#include <Eigen/Core>
using namespace Eigen;

// GeoReact includes
#include "EquilibriumCondition.h"
#include "../Multiphase.h"
#include "../ReactionSystem.h"

class EquilibriumSolver 
{
public:
	struct Options
	{ unsigned maxiterations; double tolerance; };
	
	struct Diagnostic
	{ unsigned iterations; bool converged; };
	
public:
	EquilibriumSolver();
	
	EquilibriumSolver(const Multiphase& multiphase, const ReactionSystem& eReactions);
	
	void
	Initialiaze(const Multiphase& multiphase, const ReactionSystem& eReactions);
	
	void
	SetEquilibriumConditions(const vector<EquilibriumCondition>& eConditions);
	
	void
	SetTotalMassMode();
	
	void
	SetSolverOptions(double tolerance, unsigned maxiterations);
	
	const Diagnostic
	Solve(double T, double P, VectorXd& n) const;
	
	const Diagnostic
	SolveTotalMassMode(const VectorXd& uj, double T, double P, VectorXd& n) const;
	
	const EquilibriumCondition
	FixedMass(const string& species, const FunctionTPna& h);
	
	const EquilibriumCondition
	FixedMass(const string& species, double mass);
	
	const EquilibriumCondition
	FixedTotalMass(const string& species, double totalmass);
	
	const EquilibriumCondition
	FixedActivity(const string& species, double activity);
	
	const EquilibriumCondition
	FixedAcidity(double pH);
	
	const EquilibriumCondition
	FixedPartialPressure(const string& gas, double Pg);
	
	const EquilibriumCondition
	ChargeBalanced();
	
	const VectorXd
	EstimateInitialGuess(const vector<EquilibriumCondition>& eConditions);
	
	friend class EquilibriumProblem;
	friend class EquilibriumProblemTotalMassMode;
	
private:
	/// The phase system object which is used to compute activities of the species for every phase 
	Multiphase multiphase;
	
	/// The system of equilibrium-controlled reactions that defines the set of mass-action law equations
	ReactionSystem eReactions;
	
	/// The primary species (jSpecies) and equilibrium species (eSpecies) in the multiphase system
	vector<string> jSpecies, eSpecies;
	
	/// The indexes of the primary species (jIndexes) and equilibrium species (eIndexes) in the multiphase system
	vector<Index> jIndexes, eIndexes;
	
	/// The number of primary (Nj) and equilibrium (Ne) species
	unsigned Nj, Ne;
	
	/// The stoichiometric matrix of the equilibrium-controlled reactions
	MatrixXd vej;
	
	/// The matrices that define the set of equilibrium constraints in the system 
	MatrixXd alpha, beta;
	
	/// The functions that evaluate the right-hand side vector h of the system of equilibrium-constraint equations
	vector<FunctionTPna> h;
	
	/// The tolerance for the iterative procedure
	double tolerance;
	
	/// The maximum number of iterations for the iterative procedure
	unsigned maxiterations;
	
private:
	void
	UnderRelax(VectorXd& n, const VectorXd& nOld) const;
};

#endif /* EQUILIBRIUMSOLVER_H_ */
